#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>
#include <vector>
#include "Function.h"
#include "Polynomial.h"
#include "Interpolation.h"
#include "spline.h"
#include "Curve_fitting.h"
#include <jsoncpp/json/json.h>

using namespace std;

class fun : public Function
{
public:
  double operator()(double _x)
  {
    return 1.0/(1.0 + 25.0*_x*_x);
  }
  double diff(double _x)
  {
    return -50.0*_x/pow(1 + 25.0*_x*_x, 2);
  }
  double diff2(double _x)
  {
    return 50.0*(75.0*_x*_x - 1)/pow(1.0 + 25.0*_x*_x, 3);
  }
};

int main()
{
  Json::Reader reader;
  Json::Value root;
  ofstream fout1,fout2;
  ifstream in("A.json",ios::binary);
  if (!in.is_open())
    {
      cout << "error: cannot open file." <<endl;
      exit(0);
    }
  double l, r, err;
  int _Boundary_conditions, _order;
  vector<int> Point_Num;
  fun f;

  if (reader.parse(in,root))
    {
      l = root["interval_l"].asDouble();
      r = root["interval_r"].asDouble();
      for (int i = 0; i < root["n"].size(); i++)
	{
	  int _n = root["n"][i].asInt();
	  Point_Num.push_back(_n);
	}
      _Boundary_conditions = root["Boundary_conditions"].asInt();
      _order = root["order"].asInt();
    }

  fout1.open("A_result_points.txt");
  fout2.open("A_result_error.txt");
  fout1 << "order = " << _order << ", Boundary_conditions = " << _Boundary_conditions  << endl;
  fout2 << "max interpolation error with order = " << _order << ", Boundary_conditions = " << _Boundary_conditions << endl;
  
  for (int i = 0; i < Point_Num.size(); i++)
    {
      vector<double> T;
      for (int j = 0; j < Point_Num[i]; j++)
	{
	  T.push_back(l + j * (r - l) / (Point_Num[i] - 1));
	}
      ppForm_interpolation pp( f, T, _Boundary_conditions, _order);
      pp.solve();
      fout1 << "N = " << Point_Num[i] << endl;
      for (int j = 0; j < 200; j++ )
	{
	  fout1 << pp(l + 0.01 * j) << ",";
	}
      fout1 << pp(r) <<endl;
      err = 0;
      for (int j = 0; j < Point_Num[i] - 1; j++)
	{
	  if (fabs(pp(0.5*(T[j] + T[j+1])) - f(0.5*(T[j] + T[j+1]))) > err)
	    {
	      err = fabs(pp(0.5*(T[j] + T[j+1])) - f(0.5*(T[j] + T[j+1])));
	    }
	}
      fout2 << "N = " << Point_Num[i] << ": "<< err << endl;
    }
  fout1.close();
  fout2.close();
  return 0;
}
